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ABSTRACT 


This report presents the SHORE code user's manual. The 
SHORE code computes the transient response of shells of 
revolution to initial velocities and surface pressure histories; 
it includes nonlinear material and geometric effects. 


This document presents the theoretical development of the 
SHORE code and two sample problems to illustrate the use 


of the code. 
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NOMENCLATURE 


elastic moduli 

linear strain hardening modulus 
impulse 

stress resultants 

stress couples 

transverse stress resultants 
gravitational constant 

shell wall thickness 

index in ¢-direction 

index in @-direction 

time increment index 

linear foundation springs 


surface pressures 


radius of shell in @-direction norma! to the shell wall 
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radius of shell in ¢-direction normal to the shell wall 


time 
displacement components 


radial coordinate 


mesh spacing in 0--direction at the i 


mesh spacing in the ¢-direction 


integration time step 
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displacement increments 
rotations 

total strain increments 
stress increments 

strains in the mid-surface 
circumferential coordinate 
bending and twisting strains 
damping coefficient 
plasticity proportionality factor 
Poisson's ratios 

weight density 

yield stress 

effective yield stress 


stresses 
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angle between r and shell axis (meridional coordinate) 
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INTRODUCTION 


The SHORE code computes the transient response of shells with midsurfaces defined 
by a surface of revolution. The equations solved by two-dimensional finite difference 
techniques are those derived by Sanders [1,2] and include the nonlinear geometric 
effects. Time integration is by the central difference method, Inelastic material be- 
havior based on Hill's anisotropic yield function (3] is valid for initially orthotropic 
materials for both perfect plasticity and linear strain hardening. A fracture criterion 
also is available. The loadings considered are initial radial velocity, initial step 
change in temperature (never checked out), and surface pressure histories; all of these 
can be specified for each mesh point. The shell can have either one or two different 
material layers obeying the Kirchoff assumptions with additional weight and coordinate 
direction springs at each mesh point. The boundary conditions include free and fixed 
simple supports, clamped supports, free edges, conditions of symmetry and anti- 
symmetry. and inhomogeneous formulations. The code is checked out for cylinders 
and cones for elastic, inelastic, linear, and nonlinear kinematic response. Other 
shapes have been analyzed, but no real check cases are available. 


The code development started as a extension to the work by Franke [4] and Hubka [5] 

to include more general geometry, boundary conditions, and material characteristics. 
Therefore, the preliminary versions of the SHORE code closely resemble the CYLINDER 
code [4]. As more experience was gained and more understanding developed, the 
SHORE code became independent. Now, all that exists of the original starting point is 
the framework of the overall logic and a few of the subroutine names, 


The outstanding feature of the code is the relatively fast run times compared with other 
two-dimensional shell codes. This results from using lower order difference expres- 
sions O(ds?) (As, a spatial increment) in conjunction with a checkerboard type half 
spacing grid system. This system appears to produce results equivalent in accuracy 
to those obtained by codes STAR [6] and GIRLS1 [7], using O(ds*) finite difference 


ix 
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expressions, The checkerboard grid is a slightly modified version of one employed by 
Cassell [8]. The use of the OAs”) difference expressions considerably reduces the 
number of computations per integration time step. In fact, algebraic forms are used 
throughout the coding wherever possible in order to avoid time consuming numerical 


techniques. 
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Section 1 
PHILOSOPHY 


This section contains both the philosophy of the actual code and the philosophy of the 
user's guide, These topics are probably the most important for the beginning user 
if one is to make maximum use of the SHORE code. 


The code was written to provide an efficient, relatively easy-to-use response analysis 
code for two-dimensional shell structures. This task was not as easy as it may appear 
because it has taken a few years for the SHORE code to evolve to its present form. 
This present version may not be permanent, but apparently it has more permanency 
than past versions. The greatest struggles occur in (1) finding an accurate simpie 
spatial finite differencing scheme and (2) understanding the proper implementation of 


boundary conditions. 


Spatial grids employing whole stations with O(ds”) and Olas") finite difference expres- 

sions were tried before the present half-station, OAs”) , tcchnique now used. The 

whole-station, Olas") , technique is accurate, but requires a very large amount of 

time for computation. The whole station, OAs’) , technique was not accurate when 

compared with other known solutions. For reasons not clearly understood, the half- > 
station, O(ds’) , technique has the accuracy and — because of the very simple expres- 

| sion for the derivative ~ has computational efficiency; thus, the half-station technique 

is used in the code, 


For the boundary conditions, the technique is used of evaluating fictitious points out- 
side the boundary, Although techniques involving backward and forward derivatives 
without fictitious points were tried briefly with resulting stability problems, no real 

| experiments have been conducted with the boundary condition solution technique. The 
technique of using fictitious points is so much easier, especially for the OAs”) finite 
| difference expressions, that there is no apparent contest. 
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It is believed that the user's guide is adequately detailed. Hopefully, this detail will 
enable the user to research the code and make changes to work a wider variety of 
problems than those encompassed by the production version of the code; the second 
sample problem, Section 4,2, illustrates this. The user may want to make other 
changes, ¢.g., temperature dependent properties and alternatives to the built-in 
pressure histories. 


A few remarks are appropriate concerning general use of the code. The run times can 
be quite large for many analyses; hence, it is most important that the problem is setup 
correctly. The first run should be for a very short period of time with ample printing 
of results. The input data should be carefuliy checked: is this the problem to be 
solved? The output should be checked: do the answers look reasonable? Judgment of 
reasonableness will come only through experience gained by studying the results. 

When it is certain that the problem has been set up correctly, the printed output may be 
reduced; the plots can be used for fast runs, and the analysis can be run to completion. 


Run times vary with the amount of data printed and plotted; printed data are very costly. 


The following are rough estimates for typical output. 


Run time (seconds) = (Number of meridional stations)(Number of circumfer- 


ential stations) (Number of integrations) (4.5 x 107%) for 1 layer elastic runs. 


The factor (4.5 x 10°°) becomes approximately (8 x 10°) for 4-layer elastic plastic 


analyses, For nonlinear geometry analyses, multiply the factors by approximately 1.6, 


If very short runs with a small number cf stations are made, more run time will be 
required as the setup time will begin to dominate the run time. The setup time varies 
from 10 to 40 sec on the 1108 Exec 8 system. The maximum, 40 sec, includes 
recomplilation of several subroutines. 
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Section 2 
ANALYSIS 


2.1 BASIC EQUATIONS 
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The shell equations solved by the SHORE code are those resulting from Sanders [1, 2] 
and include both the linear and nonlinear kinematics. For use in the SHORE code, the 


general equations presented by Sanders are specialized to a shell whose initial mid- 


surface is definable by a surface of revolution. This greatly reduces the number of 


derivatives in the circumferential direction without severely restricting the class of 


practical problems that can be solved. For the details, see Section 2.5 on shell 


geometry. 


The equilibrium equations are 
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The underlined terms are those included for the nonlinear 


kinematics solution. 
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(2. 1-3a) 


(2. 1-3b) 


(2, 1-3c) 


If yield occurs, the stresses are determined through an incremental plasticity theory 


which is discussed in Section 2.4. 


The stress resultants N ¢ No »N,,,M »My , and M,, —used in the equilibrium 


90’ ’ po 
Eqs. (2. 1-1a—c) — are defined as follows, 
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x °% 
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and 
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Equations (2,1-1) through (2, 1-4) constitute the complete set of equations needed to 
solve the response of a shell of revolution. The next section describes the solution 


technique employed in the SHORE code. 
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2.2 SOLUTION METHOD 


Both the time integration and spatial shell equations are solved by finite difference 
techniques. In both cases, simplicity has been emphasized and the lowest possible 
order derivatives have been used to effect a proper solution. Table 1 lists the differ- 


encing expressions used in the SHORE code. 


For the time integration, a simple procedure is employed using the central difference 
expression with displacement increments. The displacement increments are velocities, 
so instead of computing accelerations as is done in many codes, the lower order veloc- 
ities are used in the SHORE code. References 9 and 10 are excellent source material for 


for this subject. 


For a function F 


i ee aie ae 
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é6F = FL. -F, 
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r AF = F.,,- F, 
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FINITE DIFFERENCE EXPRESSIONS USED IN THE SHORE CODE 


Central 
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where 6F is the past increment in F and AF is the next increment in F. Thus, 


Eq. (2.2-1) can be written as 


oF _ AF + 6F 


Bion —2 at 2At 


Refering to the equilibrium equations (Eqs. 2.1-la—c), we can write them in the 


following form. 


a (eree _ (QF + SF) _ ph [ar - 6F\| | 
E KF a( DAt ) »( 3 )) = £00.10, 


or 


; d | ‘area: 
AF = |-6F(s% - 5) +p, -k F+ L (N,M) (a+ 4 
(A gate fof K 2at  oht 


So at any given time (k) — and with 65F (the past increment) known — the next increment 
in F, AF, can be determined. Of course, when AF is computed it becomes the 6F 
for the next increment in time, At, so one need only store the variable AF. Indeed, 
a recursion relationship for AF can be written as follows in terms of computer 


manipulation, 
AF = -ar(a - 4) +p - KF +L (NM), fee 
gdt gAt 
-— (2,2-2) 


The total value of F at anytime k is simply 


k 
F. = > OF; 
j=l 
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or expressed in a computational form 


The real advantage of using these increments in displacement will become evident 


when we deal with the plasticity theory. 


As shown by many investigators, this explicit time integration is stable for the 


At < Ax, I (2, 2-3) 


where Ax indicates the smallest spatial mesh size and Ve/eE represents the recip- 
rocal cf the maximum sound speed in shell wall. For the input data, a check is made 
in subroutine RDINPT. If input At is larger than that for stability, the code will run 
with its own value of At. In cases where additional springs or mass have been added 
to the shell, even this At may prove to be unstable, if the added terms result in a 
higher system frequency, This smaller At can be estimated from the following 


following At. 


formula, 


1/2 


1 


] 1 
p + (added weight)/thickness i p(thickness) + (added weight) 
gE 2 


gr k 


At < &x (2. 2-4) 


Experience indicates that At should be chosen as 80 percent of At ax to ensure 


stability under all conditions. 


The spatial sclution of the equilibrium equations is not as straightforward; however, 
it is simple ‘f the user has previously worked with finite differences, The following 


grid system [8] is used, 
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The main feature of this grid system is the half spacing of the u,v displacements 
between the radial, w, displacement, This provides a strong coupling between the 
three equilibrium equations not generally provided in lower order difference expres- 

sions. A problem in using this grid exists in solving the plasticity problem in | 


which it is necessary to have the stresses and strains known at the same point 


to properly apply the yield condition and compute the stress increments. This is 
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done by shifting the points for shear quantities, o, io the same pivotal points used for 
w,e@. This destroys some of the simplicity of the various required derivatives, but 


it eliminates the problem without any apparent loss of accuracy. 


To illustrate application of the finite difference formulas, the two following expressions 


are presented. 


a Wij 
€ = et fa en , (linear only) 
oa tie ee 
and 
-M M -M M 
%i-4 ae dint ; H ” cos o; $9; i-1 : 09; 41 
= - Se a ee ee We 
Ooh 2A (Tyo) ij 5 rT; sin >; 24 (r; é sin ¢;) 


evaluated for w. Continuing in this way with the mesh configuration of Fig, 2, 2-1 
and the finite difference expressions in Table 1, all of the basic equations in Section 2. 1 


can be written in finite difference form. 


By referring to Fig. 2,2-2 which shows the SHORE code flowchart, the computational 
sequence of events can be quickly understood, 


Initialization. A few quantities used in the code are initialized to avoid diag- 
nosistics relating to having a value assigned for quantities used later in the code but 
for which no value may have been computed for the case being run. No initialization 


is made that would enable one to stack cases, 


Set Up Drum. Next, the drum allocations needed for the plot routines are 


set up. 


RDINPT. This subroutine is called from SHORE, and it reads the input data 


from cards and prints out the input data. 
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INITIALIZATION 
SET UP DRUM 


RDINPT 
Read and write input data 


COMPUTE CONSTANTS 


SET INITIAL CONDITIONS 
Call Entry IDISP, Subroutine 3NDCON 


UPDATE TIME 


SET DISPLACEMENT BOUNDARY CONDITIONS 
Call ENTRY DBC, Subroutine BNDCON 


STRAINS & STRESSES 
Compute 6€ and e—- 50 anda 
Check for yielding; if yes, recompute o 


UPDATE DISPLACEMENTS 


A aaa es + Aw, etc 


STRESS RESULTANTS { 
Evaluate Equations (2.1.4) 


| 


SET RESULTANT FICTITIOUS POINTS 
Call ENTRY SRBC, Subroutine BNDCON 


NEW DISPLACEMENT INCREMENTS 
Compute Au, Av and Aw from Equations 
(2.2-2) and (2.1-1) 


Fig. 2.2-2 SHORE Code Flowchart ; 
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Constants. Various constants used throughout the program are evaluated 


from the input data. 


Initial Conditions, Either initial thermal strains or thermal excitation (never 


checked out or used yet) or the initial velocity is computed (see Section 2.6, Loading). 


Update Time. This is the beginning of the time integration loop; at this point, 


the value of k, the time index, is incremented by one. 


Displacement Boundary Conditions, The desired boundary conditions are set 
(ENTRY DBC in Subroutine BNDCON), (see Section 2.3), Because the code actually 
computes the displacement increments, the boundary conditions are set on the incre- 


ments only. 


Strains and Stresses. The increments in strain are computed based on the 
displacement increments and their derivatives. The total strain is computed from 
Ex = &-1 +5e€. The increment in stress is computed based on the elastic constitutive 
relationships, Eqs. (2.1-3). The total stress is computed from O, = 4+ 60, then 
a check on the yield condition is made. If yielding has occurred, the total strain, 
stress and the strain increments are used in the subroutine INELST, see Section 2. 4, 


to determine the proper total stress. 


In addition, a fracture criterion is available, The criterion results from Hoffman [11] 
and accounts for differences in tensile and compressive fracture stresses. If this 
criterion is exceeded, the function YIELD is set equal to -1. and the stresses at this 


point will be zero afterwards. 


Update Displacements. The present displacement increment is added to the 
sum of the past displacement increments. 


Stress Resultants. The values of the stress resultants, Eq. (2. 1-4) are com- 
puted, For elastic cases, an algebraic form of the integration through the thickness 
is used. For inelastic cases, the stresses are determined at five points through the 


thickness; Simpson's rule is used for the integrations. 
2-13 
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Plot Out. Different computed quantities are plotted at various spatial locations 
and times as specified in the input data. The subroutine PLTNG selects the desired 
data and then calls subroutine SCPLT which stores the data on a drum for actual plotting 
after the computations have been completed, After the she!l response has been com- 
puted in subroutine SHORE, control returns to BEGIN, the driving program, at which 
point PLOT! is called and the plots are prepared. The following subroutines are 
connected with the plotting: SCPLT, PLOT1, PLOT2, DXDY2, NREAD and NWRITE. 


Print Out, Various computed quantities at various times as specified in the 
input data are printed out by calling subroutine PRNT. 


Resultant Boundary Conditions. As discussed in Section 2.3, the fictitious 
points for the stress resultants are set by calling ENTRY SRBC in subroutine BNDCON. 


New Displacement Increments. Using the time integration technique outlined 
at the beginning of this section, Eq. (2. 2-2) is applied to the three equilibrium Eqs. 
(2. 1-1) expressed in finite difference form to compute the next displacement increment, 
At this point, the code returns to the Update Time statement, and the entire sequence 
is repeated to compute the next displacement increments. 


This completes the outline of the flow of computations in the SHORE code, but a few 
general comments are still in order. The production version of the code uses the grid 

with the shear quantities computed at the w pivotal points. If interest exists only in 

elastic analysis, a version of the code with the grid in Fig. 2,2-1 is available. Although 

the code is called SHORE, it is actually only a subroutine which is called from BEGIN. 

BEGIN is the main program —a very simple driver program, PLOT1 also is called 

from BEGIN. Because the code is very large there are several segments and overlays; 

these can be determined from the map obtained from a listing. Also, at each time step 

the foregoing sequence of operations is performed for each spatial point associated with 

the shell. 
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| The following nomenclature is used tn SHORE: 
SCi1 1/A(r, 9) 
i SC2 1/A6 
TATL temperature matrix 
EPP Ey 
EPT Eg 
. GPT < 40 
| SP oy 
ST Tg 
{ TTP T $9 
U,V,W u,v,w 
f DU , DV , DW Au , Av , Aw 
NP N 
p 
[ NT Ny 
NPT Ngo 
MP M 
d 
[ MT M, 
MPT Myo 
ZZ Zz 


T,TU,TV,TA,TUA,TVA thickness of inner shell layer at w,u,v, respectively 
and outer shell layer at w,u,v, respectively 


I meridional direction index 

J circumferential direction index 
K time increment index 

L thickness direction index 
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2.3 BOUNDARY CONDITIONS by Charles Rankin 


In principle, there are 16 possible boundary conditions at each end of the shell, as 


described by the following combinations which apply at s = 0 and S = L: 


Ny = 0 or u= 0 (2. 

Neo +2 (F- z Myo t BMY + Noe = 0 eee lee (2 
aM 50 

2% ‘ a(ré sin od) a ON, o ® N46 = 0 or w= Q (2 

M = 0 or ®& = 0 (2 


Additional conditions arise when using symmetry to cut down on the calculation. 
w or v symmetric, 


=0 and 2% = 9 
88 
3 
X= 9 and oy = 0 
: 0s 


Because of half-spacing, only one equation is needed for symmetric u (because 
and v are 0 at the boundary): du /ds = 0 at boundaries. 


Antisymmetric conditions for w,v, and u are 
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w =0 and — =0 


2 

v=0 and er, 
2 

0s 


SHORE can operate with the following boundary conditions at either s =0 or s = L: 


(1) Clamped (u = w = dw/ds = v = 0) 

(2) Fixed simple support (u =w = 0, M, =0, v = 0) 

(3) Free simple support (w = 0, N, =0, M. =0, v = 0) 
(4) Free edge. v = 0 (axisymmetric case) 

(5) Unrestricted free edge 

(6) All symmetry combinations of u,v, and w 


All boundary conditions, except for the two free edge cases, require that one fictitious 
point outside the shell be used to set the boundary conditions. By setting these points, 
central difference formulas expressing the derivatives in the boundary condition equa- 
tions are satisfied. Other points not required hy the boundary equations but needed for 
quantities in the equilibrium expression are extrapolated by a second-order formula, 


which for a variable y is 
Yi-n = 39 ~ 8¥in1 + Vine 


Points treated like this shall be referred to as extrapolated. 


The clamped condition requires that u = w = dw/ds = 0. Let (u, i’ v4 j’ wy, ) 
represent a fictitious point at the s = 0 boundary at gridpoint y along the circum- 


ference, and let Unc, j ; YNC+1,j : Wye +1, ;) represent the same fictitious point at 


the right boundary. (Note that there is one less u point because the u points are 
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halfway between the w points; see Fig. 2.2-1 and discussion.) Then, at s = 0 and 


s=L, 
Ww = Ww ; w Ww (or ow 0) 
Led 3,j NC+1,j NC-1,j F) 
Wj = Up j ; WNC,j WC-1,j (for u = 0) 
Wo j = 0 : YNC,j = 0 (for w = 0) 
Yo j = 0 : YNC,j = 0 (for v = 0) 


The term v, is extrapolated. 


1 


For w (or v) symmetric (w' = w™ = 0), the third derivative terms require a dif- 


ference expression shifted one point to the right (or left for S = L), and are as 


follows: 


d2w | Z (Ww, , - wg + 3wy y - Wy 4) 
3 
ds 50 As 
d°w _ _ *ne-2,j ~ *¥no-1,) * 2nc,j ~ “ne+1,j) 
3 3 
ds Seb As 


The equations for symmetric w (or v) become 


Symmetric 


7. Soa > Wneo+1,j ~ YNne-1,j 
a a 4 os ) : = ay -w ) 
= 3 4ws 5 - 4,5 ; Wne,j * 3 @¥nc-1,j ~ “No-2,j 
u requires 

Pgh 1 amet “e,j ~ “NC-1,j 
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Similar treatment for antisymmetric conditions yields for w (or v) 


we, . = 0 : WNO45 oe 


Wnc+1.j ~ ~“NC-1,j 


and for u 


Vig > les “Nc,j ~ “NC-1,j 


In the fixed simple-support case, u = 0 and w = 0 at the boundary, which requires 
that 


NC,j  —“NC-1,j 


we, . = 0 a veg °° 


In addition the requirement exists that M 6 = 0 atthe boundary, This relation forms 
a complicated nonlinear equation involving numerical integration (to find M 5 from 
stresses through the thickness) for the unknown quantity Wy ij (or ’ Wnt ij The 
regula falsi method is used to find the root. Initial estimates for Wij are found by 
requiring that the dominant term in M, (i.e., Ly + UGK o) be zero. [See Eqs. (2.1.2d) 


and (2.1.2e) for details. ] 


$ 


Only a very few iterations are required to find a good value for w 1) 
The free simple support used in SHORE is a modification of the v = 0, w = 0, 
N, = 0, and M, = 0 condition [Eqs. (2.3-1) through (2.3-4)]. Instead of setting 


p 
w = 0, the perpendicular distance from the boundary point (Ww) to the axis of the 
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shell is kept constant. Thus, if the shell translates in the axial direction (as in the case 
of a cone with free simple-support condition at each end), no changes in the shell configu- 


ration will take place. The (translationally invariant) equation for w is 


ucosd, ., = 0 


sin 
which becomes 


(u +u, .) cos ¢ 
w..es eee (ant Fp cat 4 (2,3-5a) 


for s = 0 and 


Uy, ;,7t U ) cos ¢ 
* = - -NC-1,j NC, j' NC (2. 3-5b) 


NC,j 2 sin ?nc 


at s = L. Equations N ¢ =M ¢ = 0 [together with Eqs. (2.3-5)] are sufficient to 


fix the values of Wj ' hj ; 


Here we have two nonlinear equations in two unknowns [after Eqs. (2.3-5)] have been 
absorbed] which must be solved numerically. To this end, a generalization of the 
regula falsi is used which was developed by Gauss [14] and is applicable to sets of 
equations. Although the method appears to be used infrequently and although a precise 
consideration of its convergence properties has not been carried out, the method seems 
to converge as well in many dimensions as the regula falsi (to which it reduces in one 
dimension) does in one. The advantage over the numerical versions of Newton's 
method is that only values of the equations (and not its derivatives) are needed. Coding 
for the regula falsi is in the subroutine ROOTS, which in turn needs SET to set values 
of the stress and moment resultants, STRESS to calculate these resultants, and LINEQ2 
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to effectively solve a set of linear equations for a possibly ill-conditioned matrix of 
coefficients. In case of lack of convergence, the latest values of the unknowns (usually 
near zero for the displacement increments used in SHORE) are returned and a message 


as follows is printed out: 
TOO MANY ITERATIONS 


TIME = XX MILLISEC II = XX JJ = XX NDEM = XX Kd = XX 
XXXX, XXXX. XXXK --- 


for iterations exceeding a safety value preset in BNDCON, or 


ARGUMENT VALUES OUT OF RANGE 


TIME = XX MILLISEC II = XXJJ = XX NDEM = XX Kd = XX 
XXX. XXX 


if the values for the unknown go out of a specified range (in BNDCON). 


Here XX refers to printed values of specified quantities, time is the time step when 
the error occurred, II and JJ refer to grid points (do not include fictitious points ex- 
cept for free edge case; see Section 2.3), NDEM refers to the number of equations, 
and KJ is the number of iterations. The numbers below these are the argument values 
and the function values, respectively. For asymmetric runs, a few spurious points 


should not adversely affect the results. 


For all the above (non-free-edge cases), only one fictitious point is needed. Basically, 
this results from the fact that the equilibrium equations do not need to be satisfied on 
the boundary for the clamped, simple-support, and symmetry conditions; indeed, the 
constraints imposed violate the equilibrium equations. The true solution actually 
starts one point in from the boundary. In the free case, on the other hand, the bound- 
ary is responding like the rest of the shell, and the equilibrium equation must be satis- 


fied there just as in the interior. These equations are fourth order in w cid second 
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order in u and v; hence, theoretically two fictitious points in w and one in u and 
v are needed. Thus, one additional point is added to the left (right) wherever a free 
boundary is encountered. Unlike non-free cases, one of the fictitious points is treated 
as if part of the shell: It must have a thickness, initial temperatures (if applicable), or 
any other quantity that would be required to extrapolate the shell one point beyond the 
boundary. If the free edge is on the right (I = NC), this is a simple matter of adding 
one point. If on the left, however, each input grid point for these quantities is dis- 
placed one point. [ Thus, pt.(3,J) for simple support becomes (4,J) for free edge on 
left.) Note that for R, RP, and PHI, two fictitious points are required. (See Sec- 
tion 3.) In all cases, this free edge fictitious point is treated as a real part of the 
shell with real physical input extrapolated from the interior. The plotting, on the 
other hand, is treated exactly the same, regardless of the boundary conditions. 


Two free edge choices are available, v restricted and v not restricted. The v = 0 
case is used primarily for axisymmetric jobs. If the user assigns his axisymmetric 
run to the unrestricted free edge, the program will correct the mistake. On the other 
hand, the v = 0 option is retained for asymmetric runs because of a possible cut in 
run time for cases where knowledge of v is not important. 


From here on, only the left boundary will be discussed; the right is treated similarly. 


When v = 0, the points in question are Vo j (I = 3 at the boundary), Yj" Wo 5’ 


ha Ee pte on ey 


The terms V4, j and u, J contribute but little (just as in the case of Ya, j for clamped 

case) and are extrapolated; Vo J is set equal to -v 4,j (which in any case is zero for 
axisymmetric runs), thus treating v as "antisymmetric"; Wo | i’ Wa j’ Up j form ‘ 
three unknowns for the three equations 


N, =0 (2.3-6a) 
M, = 0 (2.3-6b) 


OM. 6 


Qy +3 sin }) ~ ByNyg = 0 (2.3-6c) 
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at the boundary. Because Eqs. (2.3-6a) and (2.3-6b) do not contain Wy \’ they can 
be solved for the pair of unknowns Wo i’ Uy j using ROOTS. Then, Wy j is solved 


using Eq. (2.3-6c) after Wo j and Us j are known. 


In all cases, initial estimates for the unknowns Wo j and Uy j are taken from 
dominant terms in M e and N a respectively: 

Wo ij from Ky + V 6k = 0 

Wo ij from €% + Voeg = 0 


The guess for Ww) j is extrapolated. 


The unrestricted free edge is treated similarly, except one more equation comes 


into play: 
1/3 1 1 ey = 
Ngo 3 (7 Myat 2 Met = 0 (2.3-6d) 


Equation (2.3-6d) is included with Eq. (2.3-6a) and (2.3-6b) for the unknowns Uy i’ 


and w and v are extrapolated. Because cross terms in 


Veg + MAG 
2,j Rs lana | 1,j 
Eqs. (2.3-6) weakly couple all terms around the circumference, these three equations 


must be solved more than once around the circumference. They are solved for each J 
value one at a time around and around the shell until convergence is attained. In prac- 


tice, only two rounds of the shell are needed; if the number of circuits exceeds three, 


the program prints out 
OUCH, IT HURTS, KITER IS GREATER THAN 3. HELP, HELP 


and proceeds merrily along its way to solve for Ww) j from Eq. (2.3-6c) as in the 


restricted free case. 
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In the circumferential direction, a plane of symmetry exigis at 6 = 0 and O° =f) 
that the displacement must be set as follows for all stations I along the length of the 
shell: 


e Symmetric loading 


Ca ; 4 NBt1 ~ “4, NB-1 
Nig 2 ie , Yi NB “i, NB-1 
hs te Wa ; Wi NBt1 ~ “i,NB-1 


e Antisymmetric loading — Reverse signs on right-hand side. 


The displacement boundary conditions are set in the BNDCON subroutine at the entry 
point DBC. Because the code actually computes displacement increments in u, v, 
and w called DU, DV, and DW in the code, these variables are seen in ENTRY DBC. 
CIRC computes the displacement boundary values for 6 = 0 and @ = a 


Entry SRBC is called after the stress resultants are calculated to set M 6 and N é 

at zero wherever applicable. Fictitious values of the resultants also are calculated, 
but these values are needed only when the equilibrium equations are to be solved for 
Wo jNe, ;) and Vo, i'Nc, ;) , which include only the symmetry options. All values 
are calculated by extrapolation, except for the case symmetric v, w and antisym- 


metric u, whence Eqs. (2.1-2) also are symmetric or antisymmetric, resulting in 


“Oi ” “bs j One nn One) | 
"61,5 : me 3 nots ‘ a eer 
rer : "oy “eee, j ” "ONG 
aie a BCA es ~ ne 


‘ 
[cape al —— reed 
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N = -M 
9 9no+1,j 


M = -M 


M oie 
9 ?no+1,j %9No-1,; 


| J ow 
For antisymmetric w.v., and symmetric u, reverse the signs on the right-hand side. 
SRBC calls entry RCIRC which sets fictitious stress resultants according to the above 


symmetry condition for all points along the shell at 6 = 0 and On ° 


Inhomogeneous boundary condit‘uns can also be treated with SHORE. For the axisym- 
metric case. input the axisymmetric free edge boundary condition; for asymmetric 
runs, input the free edge. Two options are available. Either the user can alter sub- 
routine SET (see next paragraph) or springs can be attached to the boundary. Springs 
treated in the boundary routine calculate an externally applied force proportional to the 
displacement in the direction of the spring, opposing motion in that direction. 


These springs are inputted in the arrays KP(I,J) (meridional springs), KC(?,J) (circum- 
ferential springs), and KZ(I.J) (radial springs). No bending springs are available. 
Springs on the left occupy KP(i,J), KC(1,J), and KZ(1,J) and on the right occupy 


KP(NDP+2+ L, J) 
KC(NDP+2+L,J) 
KZ(NDP+2+L, J) 


where J refers to any position on the circumference, NDP is the number of meridional 
shell segments, and L = 1 if there is also free edge on the left, and zero otherwise. 


If one wants to put an arbitrary force on the edge, he calculates these forces and makes 
sure their values are available to SET at statement 1251 (sequence No. 48), and then 
subtracts these values somewhere in the statement beginning with seq. No. 48 (meridi- 
onal), No. 50 (bending), and No. 52 (circumferential shear). Radial shear information 
must be made available separately at or after statement No. 1256 (seq. No. 58), and 
the statement from which the force must be subtracted begins at seq. No. 59. 
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Note that I3 = 1 onthe left boundary and I3 = -1 on the right boundary; this indi- 


cates the side that one is on when both edges are free. J refers to 


the point along the 


circumference, and COMMON/PIG/PI,G, TIME makes available the integration vari- 
able TIME (already available to SET) to user-supplied subprograms. TIME is in 


seconds. 


The following example of change cards for SET is intended to show 
many ways a user can input time-dependent forces on boundaries. 
that a function is supplied by the user as follows: 


FORCE (13,J,L) = force or moment at edge 
I3 < 0 right boundary 
I3 > 0 left boundary 
J- point on circumference 
L=1 meridional force 
L = 2 bending moment 
L=3 circumferential shear 
L=4 radial shear 


The changes in SET become 


@FOR,XS SET,SET,SET 


-49 

: *-FORCE(I3, J, 1) 
-50 

*-FORCE(I3,J, 2) 
-56 

*-FORCE(I3, J,3) 
-63 


*-FORCE(I3, J, 4) 
@MAP,S AMEN, AMENA 
@XQT AMENA 

user's data 
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2.4 PLASTICITY 


Theoretical Development. 


Two plasticity formulations are available to users of the SHORE code. The first is 
based on perfect plasticity, and the second incorporates linear strain hardeniug. Be- 
cause the code accomodates orthotropic elastic properties, the plasticity theory is 
based on Hill's anisotropic yield function, (Everything in this section is based on 
Ref. [3].) For the two-dimensional shell under consideration, Hill's yield function 


becomes (a form proposed in Ref. [2]) 


2 2 2 
o C09 Cy  F 
2F = ~-2%, 3, = 1 (2, 4-1) 


where 
o, = yield stress in the -direction 
Tn = yield stress in the 0-direction 
C12 = yield stress in shear 


Hence, anytime 2F exceeds 1, there exists a yielding for the perfectly plastic material; 
for the linear strain hardening case, Eq. (2.4-1) is modified slightly. With the above 
definition of the yield surface, it is possible to derive the proper constitutive laws for 


the plasticity theories under consideration, 
At any point in the material, the total strain can be written as the sum of the elastic 
strain and the plastic strain. Because the code computes increments in strain, the 


increments used will be in total, elastic and plastic strain, hence. 


be See Hae (2. 4-2) 
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The code does not make this incremental formulation necessary; on the contrary, the 
plasticity theories are based on increments in plastic strain. The increment in total 
strain is known from the satisfications of the equilibrium equations, and it is desired 
to determine the increments in stress, 


For plasticity, it is known that the plastic strain components are normal to the yield 
surface and proportional to the stress increments. This is expressed as 


$2 ee (2, 4-3) 
ij 80; 


for the perfectly plastic case, where y incorporates the proportionality and stress 
increments. Now the increments in total strain can be written as 


r 


Oey = E, - Vp 505) + yp ue + 604) - za, 0 + 505) 
1 


co 


ter aee 5 2 sees 2 
b€g = E, (505 V4 604) + uD (, + 605) ae oy + 60,4) (2. 4-4) 
T5 1°2 
ees a 
8¥49 = 36° * U2 (Tee * OTy9) 
"12 
There are four unknowns —éa $ 50, ,OT $0" and \ ~—s0 a fourth equation is needed; 


this is the yield function. 
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2 2 2 
o4 Oo, yy oes. : (or alee (90: ae 
| NE Se a Bis %%12 


(2. 4-5) 


Equations (2. 4-4) and (2.4-5) constitute a set of nonlinear algebraic equations; to solve 
them exactly, an iteration technique is used. From physical considerations, it is 
known that the correct » is the smallest nonnegative positive number satisfying the 
above set of equations. Consequently, the iterative technique progresses as follows, 
First, guess p = 0, then solve Eq. (2.4-4) for 60, 1 50p 574 4 , and then determine 
the sign of Eq. (2.4-5), Next, guess p = 1, solve Eq. (2.4-4), and then determine 
the sign of Eq. (2.4-5), If a sign change has occurred, then 0 < p< 1; if not, the next 
choice for p is 2 times the last ». When a sign change occurs, an interpolation 


scheme is used to find the exact value of yp and 6c,, 50% : 874 9 ‘ 


For most cases of practical interest, an approximation can be made in Eqs. (2. 4-4) 
and (2, 4-5) to linearize the equations. Then, using Cramer's rule, it is possible to 
solve directly for 60 ? 50, , OT 60° and yz. The only time this technique will produce 
a completely erroneous result is if yielding occurs on the very first time step. This 
is easily noticed because a zero stress is computed for the printout of the very first 
time step. To linearize the equations, all terms quadratic in stress increments are 
dropped. Then, the following matrix equation is obtained. 
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v 20 o 
6 
“EO 0 —2 cae 2 oa, 
ts) oy 1-2 
1 205). 966 
E, 0 ey c 0,0, 50% 
6 0 lg 
2 
1 0 
0 Tad oT 
2G O49 o6 
20 (oj 2r 
“a. _%9@\ 790 0 : 
2 T1% 2 
2 %%12 
t 
be 
d 
t 
5€5 
t 
5€ 49 
ae oe o,0 a 
1-6 6 , 6%) _ oo 
on on 01% 2 
1 2 i 
(2. 4-6) 
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For the linear strain hardening theory, a very similar procedure is used, but it is 
slightly simpler in that » is now known, and the stress increments are the only un- 
knowns. For the linear strain hardening, we use the same yield function Eq. (2, 4-1) 


and the strain hardening is assumed to be of the form 


(2. 4-7) 


where % is the current value of the effective yield stress and H is the strain harden- 
ing modulus. (Note, the initial value of is equal to the effective yield stress based 
on the von Mises' condition.) After first yielding the 1 in Eq. (2.4-1) is no longer used; 
instead, 2F must exceed its last maximum before yielding occurs. By manipulating 


several equations in Ref, (3] it can be determined that 


2 


Co 
w= gay 62F) 
where Y = current value of 2F 
hence 
o 
p __o /{&(2F) 
bei; = GHY (32 6(2F) (2. 4-8) 


In this formulation, quadratic terms in the stress increments are dropped and no exact 
solution has been programmed. Before each solution, the portion of strain increment 
required to bring the stresses up to the yield surface is subtracted, so only the plastic 
increment of strain is used. [Note Eq. (2.4-8).] The following is the resulting set of 


linear equations for the stress increments 
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2 2 2 2 2 2 
0 o 0,0 oO v o. [50,0 o Co 
pe Gad aa ae a, ||P) are eit a Be ee eg | Fe 
E YH 4 3 22 Ci) E YH 2.2 3 3 7] 
fo) 0; 91% 40155 o 40) 05 20105 201% 
ge o,.T 0, T 
Pa ae 99 __ 6 g6 OT = 6¢ 
4 Bee eae 2 60 o 
Fig “9% 9% 16 
v oe 50,0 on oe a a 0,0 an 
Seo aad LY ON eYos Fa) (ge ee 5. Ot ae éc 
Ey YH 4 2 2 2 3 9 3 d Ey H 4 3 4 22 6 
Ola = Fete o> TM NE es: 5 
on OnT 1 ej 
leo ( 2098 . Soto at 
YH 2 2 9 2 6 6 
%%19 01% 9019 
(2. 4-9) 
and 
2 
re) as 70 pea) 70 
= oe, 60 
YH \ 2 2 sere 9 2 6 
79°12 7150019 vale 01%" 10 
ae 
fe Sie aes Mi ge t 
*13G * wr’ 4 |9Toe = €ge 


Computation Subroutine INELST) 


The options in INELST are INPB = 1, 2, or 3. Option 1 corresponds to a solution of 


Eq. (2.4-9). Option 2 corresponds to a solution of Eq. (2.4-6). Option 3 corresponds 


to a solution of Eqs. (2.4-4) and (2.4-5) The following nomenclature is used in INELST. 
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DEPP =! be, 
DEPT = — bet 
DGPT = 6 “ 6 
SSX = 0 ¢ 
SST = %% 
TXT = T PY, 
SZ2 = o 

HB = H 
All, Al2, etc. are the elements of the matrices 
REXB = 1/E m 
RETB = 1/ Ey 
REXNUB = Vp/ E rf 
RETNUB) = v 6! E, 
S1B = o) 
S2B = % 
$12B = V9 
MU = Ll 
YIELD = Y 
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2.5 SHELL GEOMETRY 


The basic shell element with the coordinate system and the stress resultants is shown 


below. 


and 
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The overall reference surface geometry of the shell is illustrated below 


where 


s = arc length along the meridian, 
r = radius of the shell in the @-direction measured normal to the shell 


midsurface, 


@ = angle between the axis of the shell and r, 


r $ = radius of the shell in the ¢-direction measured normal to the shell 
midsurface, 
6 = angle defining the circumferential location 


The shell wall thickness is defined as distance from the reference surface as follows. 
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For the one layer shell: 
SURFACE 


REFERENCE = — 


Z, = NEGATIVE 
Z, = POSITIVE 


For the two layer shell: 


= NEGATIVE 
= POSITIVE 
= POSITIVE 
= POSITIVE 
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2,6 LOADING 


Two loading conditions are available to the users of the SI{(ORE code. One is an im- 
pulse or initial velocity condition; the other is a radially di:ected surface pressure 
loading which, in general, represents a blast loading on a shell, 


The impulse condition is based on an impulse/unit area resulting in an initial radial 


velocity, gw/at. This is implemented in the code by prescri!ving the value of Aw 


at the end of the first time step, At. From 
_ dw 
Av = At 
and 
I = mAv 
is obtained 
_ lt 
Aw = =a 


The value of I is specified in the input data, and the Aw is evaluated in the BNDCON 


subroutine at the IDISP entry point. The formulation is 


where the minus is used because a positive impulse is in the negative w direction. 


The pressure loading is based on an approximation to a planar blast wave traversing 


the shell as some incident angle, 7, (see illustration), 
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VELOCITY OF BLAST WAVE WITH 
RESPECT TO THE SHELL VELOCITY 


The formula for computing the radial pressure P, is, 


se iy p_ (1 + cos 6) + p, (1 - cos 6)) expl-(- time - intercept time) (2. 6-1) 
Zz 2 W; A 7; 


p.. = peak windward pressure 
Py = peak leeward pressure 


t = time value used to specify decay 


The pressure loading is computed in the subroutine SURLDS, For points behind the 
wave front, the pressure is com) uted accurding to Eq. (2.6-1) and for points in front 
of the wave front, the the surface pressures are zero. The following notation is used 
in SURLDS: 


PSI = y 


THT @ angle describing the circumferential location 
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X1 == _distance from the initial contact point of the blast wave to the point on the 
shell at which P, is being calculated 


TJ = intercept time 


VBT = distance the blast wave has travelled over the shell 


PZT = P, 
PW = Py 
PL = Py 
TAU = 7+ 


Note that Py Pps and + can be arbitrary along a meridian, but the circumferentia’ 
variation is controlled by Eq. (2.6-1). 
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Section 3 


INPUT DESCRIPTION 


3.1 INPUT DATA INSTRUCTIONS 


Cards 1 and 2 (12A6) Heading Information 


(The first heading card is printed on the plots.) 


Card 3 (1215) 


Vv 


LMSC-D244589 


NDP — number of segments in the ¢-direction = 30 - NFREE 


NDT — number of segments in the 6-direction < 18 


NLAYS -— Number of segments through the thickness 


= 1 for one layer elasticity 
= 4 for one layered shell (NSL=1), multilayered elasticity 
= 5 for two layered shell (NSL=2), multilayered elasticity 


= 


IA 1A IA IA 4A OIA 


lA 


6 


oo 8 83 6 HhUOOUCOSD 


IA IA tA WA IA IA IA 


IA 


KIMP — Loading Code 
< 0 surface pressures from subroutine SURLDS 
= 0 arbitrary impulse loading for each grid point 
= 11,6) = Is. sin (78/L) cos 0 0° 
=0 90° 
= 21,8) = I, sin (78/L) 0° 
=0 90° 
= 31@,98) = I, cos @ 0° 
= 0 90° 
= 41,6) = I, 0° 
= 0 90° 
ITEMP' -— initial temperature distribution 


Vv 


0 no distribution 
0 input distribution 
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IKP — springs in ¢-direction 
= 0 no springs 
> Oirput springs 
IKT — springs in 6-direction 
= 0 no springs 
> 0 input springs 
IKZ - springs in radial direction 
= 0 no springs 
> 0 input springs 
IMS - extra weight distribution 


= 0 no distribution 

> 0 input distribution 
IAX — symmetry indicator 

= 0 asymmetric 

= -1 ring code (NDP = 0) 

= 1 axisymmetric (NDT = 0) 
NSL —- shell layer indicator 

= 1 one shell layer 

= 2 two shell layers 


NON — Kinematics indicator 
= 0 linear kinematics 


= 1 nonlinear kinematics 


In this and all subsequent cards, NFREE = number of free edges (0, 1, or 2) in 


problem. 


Card 4 (1215) Boundary Condition — KBC (10) 


KBC (1) through KBC (4) apply to the boundary at the left side of shell, S = 0. KBC(5) 
through KBC (8) are exactly like KBC (1) through KBC (4), except that they apply to 

the right side of shell, S = L. KBC (9) and KBC (10) apply to @ = 0° and 6 = THD, 
respectively. 
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The only permissible values for KBC (1) through KBC (8) are shown in the following 
table. 


KBC's Boundary conditions to be specified 
1234 (S symmetric and A antisymmetric in u, v,w) 


1111. clamped edge 

13 11 simple support (no u displacement) 

13 1 2 simple support (no edge meridional force) 

1313. freeedge, v constrained (used for axisymmetric runs) 
1314 _ freeedge, v free 

2113 Sw,v A u (use when case to be solved has plane of symmetry) 
1221 Aw,v Su 

2121 Sw,u Av 

2111 Sw Au,v 

2123 Su,v,w 

1211 Au,v,w 

1213 Aw,uSv 

1223 Aw Su,v 


For KBC(9) and KBC(10) 


1 symmetric u,w and antisymmetric v 
= 2 antisymmetric u,w and symmetric v 


Card 5 (1215) Print out flags 


IU ~ meridional displacements 
IV - circumferential displacements 
ISX — meridional stress 
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IST - circumferential stress 
ITTX ~ shear stress 

IEPX — meridional strain 

IEPT — circumferential strain 
IGXT — shear strain 

IPZT — radial pressure 

IPPT — meridional pressure 
IPTT ~ circumferential pressure 


If = 0 no print out of indicated quantity 
# 0 print out of indicated quantity 


Card 6 (8E10.3) 


x — length of shell in @-direction (arc length) (in.) 
RHO _—— density of shell wall (1b/in. °) 

NUP - ¢-direction Poisson's ratio 

NUT — 6-direction Poisson's ratio * 


| 

eo) 

oO 
| 


angle of terminal boundary (deg) 
LMD — damping parameter (b-sec/in. °) 


Card 7 (8E10. 3) 


EX —- elastic modulus in -direction (psi) 
ET — elastic modulus in 6-direction (psi) 
GG — shear modulus (psi) 

81 —- yield stress in @-direction (psi) 

$2 - yield stress in 6-direction (psi) 
S12 — yield stress in shear (psi) 


ALPHAX -— coefficient of thermal expansion in ¢-direction 
ALPHAT -— coefficient of thermal expansion in @-direction 
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Card 8 (2E10.3,15) 


SZ — effective yield stress (psi) 

H — strain hardening slope (psi) 

INP - 1 for strain hardening, 2 for perfect plasticity and 3 for ''exact" 
perfect plasticity 


Card 9 (3810. 3) (all positive numbers) 


SIC — compressive fracture stress in ¢-direction (psi) 
S1iT — tensile fracture stress in ¢-direction (psi) 
S2C — compressive fracture stress in @-direction (psi) 
S2T — tensile fracture stress in 6-direction (psi) 


$128 


fracture stress in shear (psi) 


Card 10 (8E10. 3) Only used for 2 layer shell 


EXA — elastic modulus in ¢-direction (psi) 
ETA — elastic modulus in 6-direction (psi) 
GGA — shear modulus (psi) 
S1A — yield stress in ¢-direction (psi) ones ever 
S2A — yield stress in @-direction (psi) 
S12A — yield stress in shear (psi) 
Card 11 (2E10, 3,15) Only used for 2 layer shell > 
SZA — effective yield stress (psi) 
HA — strain hardening slope (psi) 
INPA — 1 for strain hardening, outer layer 


2 for perfect plasticity 
3 for "exact" perfect plasticity 


+ 
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Card 12 (8E10. 


S1AC 
S1AT 
S2AC 
S2AT 
S12AS 


Card 13 (8E10. 


RHOA 
TA 

NUTA 
NUPA 


Card 14 (8E10. 


IMP 
VB 
PR 


Card15 (8E10. 


DT 

i 
PRINT 
ESPT 
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3) Only used for 2 layer shell 


compressive fracture stress in ¢-direction (psi) 
- tensile fracture stress in ¢-direction (psi) 


~ compressive fracture stress in 6-direction (psi) outer layer 


tensile fracture stress in 6-direction (psi) 


- fractire stress in shear (psi) 


3) Only used for 2 layer shell 


— density of shell wall (b/in,. ) 
— thickness of shell wall (in. ) 


, : : outer layer 
- o-direction Poisson's ratio 


0-direction Poisson's ratio 


I, for impulse distribution (psi-sec) 
— relative velocity of blast wave front to shell (ft/sec) 


angle between wave front and shell axis (radians) 


eo 
— 


—- integration time step (sec) 
total time (sec) 
— printing time intervals (sec) 


— time left in run before plotting will start (sec) 


Card 16 (E10. 3,15) 


PLT 
NPLT 


plotting time intervals (sec) 


0 no plots 


iH 


5 miniature plots only (3 per frame) 


7 full and miniature plots 


8 full frame plots only 
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Card 17 (8E10. 3) 


R(1) — (radius — the distance from the axis normal to the shell wall) (in. ) 
(1) If positive R is constant 


(2) If nonpositive R is a variable and each value is input (see below) 


RP (1) — radius of shell in ¢-direction (i.e. for a cylinder or cone use a large 
number like 10!) Gn. ) 
(1) If positive RP is constant 


(2) If nonpositive RP is a variable and each value is input (see below) 


T(1) — thickness of shell (in. ) 
(1) If positive T is constant 
(2) If = 0.0, T is variable in ¢-direction only and the Z array must 
be read (see below) 
(3) If <0, Tis variable in ¢ and @ direction and the Z matrix must 


be read (see below) 


PHI(1) | — angle between R and shell axis (i.e. PHI = 1/2 for a cylinder) (radians) 
(1) If positive PHI is constant 
(2) If nonpositive PHI is a variable and each value is input (see below) 


Cards 18 through 27 require special attention if a free edge is used. The shell has one 
fictitious point (or two, if two free edges are used) outside the shell which must be 
treated in every respect as if they were part of the shell itself. This includes tempera- 
ture distribution, for example. Thus, if a free edge on the left is called for, TXTL 
(L,I,J) with values starting with I = 1 at the left shell edge becomes TXTL (1, I+: 1, J), 
and TXTL (L,1,J) becomes the fictitious value needed for the boundary condition routine. 
(See Section 2.3.) 


In addition, Cards 18, 19, and 21 need one more fictitious value. For all but the free 
end conditions, this means that one value beyond the boundary of the shell is needed for 


the boundary value ro‘itie. For free end conditions, two fictitious points beyond each 


free boundary are needed. 
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Card 18 (8E10.3) Reads NDP+3 values of R if R(1) < 0.0 


Please note that the first and last radii will be past the ends of the shell. You must 
input a fictitious point beyond the ends of the shell. Just extend the shell and compute 
the R as if the shell were really there. These values are needed for the boundary 


condition computations. 


Card 19 (8E10.3) Reads NDP+3 values of RP if RP(1) < 0.0. 


Please note as in Card 18 a fictitious point is required at each end of the shell. 


Card 20 (8E10.3) (Thickness of layers) 
If T(1) = 0. 


NDP+1+NFREE values of Z(1,I,1), 
NDP+1+NFREE values of Z(2,I,1), 
NDP+1+NFREE values of Z(3,I,1), 
NDP+1+NFREE values of Z(4,1,1) are read 


only for the 2-layer shell 


If T(1) < 0.0 (Matrix Format) 
NDP+1+NFREE, I values and NDT+1, J values of Z(1,I,J) and Z(2,I,J) and Z(2,I,J) 
(plus Z(3,1,d) and Z(4,I,J) for the 2 layer shell) are read with the matrix read routine. 


I is the ¢-index 
J is the 0-index 


(Note: fictitious points needed for free edges only.) 


Card 21 (8E10.3) Reads NDP+3 values of PHI if PHI(1) < 0.9 (radians) 


(All fictitious values required. ) 
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Card 22 If KIMP < 0 
reads (8E10.3) 


(PW(I), I = 1, NDP+1+NFREE) 
(PL(I), I = 1, NDP+1+NFREE) 
(TAU(I,) I = 1, NDP+1+NFREE) 


If KIMP = 0 (Matrix Format) 
reads IXT (I,J) — psi-sec for each mesh point 


(Free edge fictitious values only.) 


Card 23 (Matrix Format) If ITEMP > 0 


read TXTL (L,I,J) — initial temperatures for each point including thickness. (L = 1 
inner fiber) 
(Free edge fictitious values only.) 


Cards 24 (Matrix Format) If IKP > 0 


reads KP(I,J) — meridional springs (1b/in.°) 

(Left and right free edge fictitious values can be used to input spring boundary condi- 
tions. Input 0 for left and right free edge fictitious values if no boundary springs are 
desired. No fictitious values needed for non-free boundaries. ) 


Cards 25 (Matrix Format) If IKT > 0 


reads KT(I,J) — circumferential springs (Ib/in.°) 
(Same as Cards 24.) 


Cards 26 (Matrix Format) If IKZ > 0 


reads KZ(I,J) — radial springs (Ib/in. °) 


(Same as Cards 24.) 


3-9 


LOCKHEED MISSILES & SPACE COMPANY 


A ee —— So - 


LMSC-D244589 


Cards 27 (Matrix Format) If IMS > 0 


reads XMS(I,J) — extra weight (Ib/in.?) 
(Free edge fictitious values only; input data required, but the values are not used in 


the computations.) 


Cards 28, 29, and 30 ar? not read if NPLT = 0. 


Card 28 (1215) 


NPW 
NPU 
NPV 
NPEMO 
NPEMI 
NPECO 
NPECI 
NPETO 
NPETI 
NPSMO 
NPSMI 
NPSCO 
NPSCI 
NPSTO 
NPSTI 
NPMM 
NPMC 
NPMT 
NPNM 
NPNC 
NPNT 


Plotting Flags 


~ Number of radial displacements to be plotted 


~ Number of meridional displacements to be plotted 


Number of circumferential displacements to be plotted 
Number of outer fiber meridional strains to be plotted 
~ Number of inner fiber meridional strains to be plotted 


— Number of outer fiber circumferential strains to be plotted 

— Number of inner fiber circumferential strains to be plotted 

— Number of outer fiber shear strains to be plotted 

-— Number of inner fiber shear strains to be plotted 

— Number of outer fiber meridional stresses to be plotted 

— Number of inner fiber meridional stresses to be plotted 

— Number of outer fiber circumferential stresses to be plotted 

— Number of inner fiber circumferential stresses to be plotted 

— Number of outer fiber shear stresses to be plotted 

— Number of inner fiber shear stresses to be plotted | 
— Number of meridional moment resultants to be plotted 

— Number of circumferential moment resultants to be plotted 
— Number of twisting moment resultants to be plotted 

— Number of meridional stress resultants to be plotted 

— Number of circumferential stress resultants to be plotted 
— Number of shear stress resultants to be plotted 


The sum of the flags cannot exceed 300. That is, only 300 unique plots can be made. 
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Cards 29 (1215) NL-Plot Array 


This reads in the mesh points at which the variables indicated on cards 28 are to be 
plotted. 


For example, say NPW = 2 and all others are zero. 


NL(1) = @ mesh location at which the first w is to be plotted 


NL(2} = 0 mesh location at which the first w is to be plotted 
NL(3) = ¢ mesh location at which the second w is to be plotted 
NL(4) = @ mesh location at which the second w is to be plotted 


3.2 MATRIX FORMAT 
Data Card Format (313, 6E10. 3) 


NR — row number 
NC — column number 
NL — layer number 


NR,NC,NL 

NR, NC+1,NL 

NR, NC+2,NL 6 values in row NR, layer NL, and 
NR, NC+3, NL starting with column NC 

NR, NC+4,NL 

NR, NC+5,NL 


> > > > PP > 


For a two-dimensional array NL = 1 


A BLANK CARD OR ZERO ROW NUMBER 
TERMINATES INPUT OF THE MATRIX 
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3.3 DIMENSION CHANGES IN SHORE 


The standard dimensions of the code are set up for a maximum of 


NDP = 30 
NDT = 18 
NLAYS = 5 


If a ring problem or axisymmetric loading case with a finer mesh is desired or some 
other arrangement of dimensions, this can be changed easily. The following parame- 


ters are to be used: 


NA = NDP + NFREE + 1 

ND = NDT + 1 

NF = NA + 2 

NG = NDT + 3 

NH = NLAYS + 1 (NLAYS = 1, 4, and 5 only), 


where NFREE = number of free edges in the problem. 


The following control cards for the 1108 EXEC 8 are needed: 


@PDP, LF,, TRY TRY 


-66 ,66 


b’ 


PARAMETER NA = XX, ND = XX, NF = XX, NG = XX NH =X 


@ADD, P CHNG 


Element CHNG, which resides on an ELT symbolic element in the SHORE program 
packet, contains the following control cards which are executed when the ADD, P CHNG 


is encountered: 


@FOR BNDCON, BNDCON, BNDCON 
@FOR CIRC, CIRC, CIRC 

@FOR SET,SET ,SET 

@FOR STRESS, STRESS ,STRESS 
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@FOR SAVE1,SAVE1,SAVE1 
@FOR INELST , INELST , INELST 
@FOR PLTNG. PLTNG, PLTNG 
@FOR PRNT,PRNT,PRNT 
@FOR RDINPT , RDINPT, RDINPT 
@FOR SHORE, SHORE ,SHORE 
@FOR SURLDS,SURLDS,SURLDS 


Each of the programs in this list is recompiled to change the dimensions. 
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Section 4 
SAMPLE PROBLEMS 


4.1 AXISYMMETRIC ELASTIC-PLASTIC CYLINDER 


This is a relatively simple problem that illustrates the strain hardening plasticity 
option and the fracture criterion. The geometry material properties, and modeling 
of the shell are described below. 


Shell Geometry: 


radius = 3 in. 

length = 6 in. 

simple support boundary conditions loading, axisymmetric and uniform = 0. 058 
psi-sec (impulse) 


Material Properties 
modulus = 10. x 10° psi 
yield = 40. x 10° psi 
strain hardening modulus = 1. x 10° psi 


fracture stress = 60. x 10° psi 


Modeling: 


20 segments along the length, time step = 1 x 10° sec 


The first listing illustrates the data cards; the second listing illustrates the printed 
output followed by the plotted output. 
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4,2 ASYMMETRIC CYLINDER (NONLINEAR ELASTIC-PERFECTLY PLASTIC) 


This sample problem is slightly more complex then the proceeding one because the 
loading is asymmetric and the code is slightly modified to take advantage of symmetry. 
The object of this computation was to compute the permanent deformation of a cylinder 
tested by SRI [12!. The loading is uniform along the length so the problem is symmetric 
about the midplane. ‘Therefore, the shell length is input as half the actual shell length, 
and the boundary conditions are specified as symmetric at the midplane. In addition, 
the subroutine BNDCON is modified slightly so that the stress resultants also obey 

the symmetry conditions at the midplane. This is seen in the listing of the input data. 
The modification to BNDCON is not absolutely necessary, but from experience it has 
been found that this produces a more accurate satisfication of the boundary conditions. 
This illustrates the philosophy discussed in Section 1 that one must study the answers 
obtained and be prepared to ''nudge" the code into the best solution obtainable. 


The shell geometry, material properties, and modeling are listed below. 


Shell Geometry: 


radius = 3 in, 

length = 6in. (actual) —3 in. (input) 

clamped boundaries at x = 0 

symmetry at x = L/2 

loading impulse = 0,122 cos 6 ~ psi-sec (0=6= 1/2) 
= 0 (r/2<6<n) 


Material Properties (Elastic-Perfectly Plastic) 


modulus = 10. x 10° psi 


shear modulus = 4, x 10° psi 
uniaxial yield = 45, x 10° psi 
shear yield = 26, x 10° psi 
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Modeling 


10 segments along the length 
18 segments along one half the circumference 
time step = 1, x 107° sec 


The first listing illustrates the input data cards. The second listing illustrates the 
printed output (not all the print out times are shown) followed by the plotted output. 
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INTERDEPARTMENTAL COMMUNICATION 


10 SHORE Users Dem) BOG PANT/ pare 2/19/73 
Phili lerwood 52-3 : Lise 

FROM iLlip Une pert 52-33 gin. 205 PLANT/ 2 en 45521 

SUBJECT: Supplement I - User's Guide to the SHORE Code [1] 


Tntroduction 


This supplement to reference [1] is to provide additional information to users 
or the SHORE code with regard to recent chanyes and added capability. The changes 
and additions are: 1) more points through the thickness for two-layered shells in the 
plastic re-sixe, 2) three-dimensional snapshot plots of all computed response quentities, 
and 3) temperature-history forcing functions with material properties a function of 


temperature. 


First there is a brief discussion of the new features and then the changes in the 
input data are presented. Items 1 and 2 are available on the file PSU*SH@REG. and 
items 1, 2 and 3 are all available on the file FGU*SH@REF « Piease note ail old input 


decks will work with the new versions of the code. 


Discussion 


The new formulation for two layered shells was implemented because many shells of 
interest have a very thick outer layer as compared with the inner layer. So now for 
the inner layer the stresses are still computed at three equally spaced stations (two 
sublayers) through the thickness. For the outer layer the stresses are now computed 
at five equally spaced stations (four sublayers) through the thickness. Hopefully, 
this increase in the number of stations throuzh the thickness will increase the accur- 
acy of the integrations performed to determine the resultants for cases involving two 


layered shells. 


The three-dimensional plotting capability is provided by an LMSC developed package, 
DRAW3D [2]. The set up for using DRAW3D is done in subroutine THREEP which is called - 
at tne times tne snapshots ere desired. The capabilities of DRAW3D ere quite versatile, 
but in THREEF the calling statement for DRAW3D has been specified for only one viewing 
angle. Aliso only the radial displacement plot has the accompanying minimum and maximum 


cross-sections. Anyway if the user wants to change the three-dimensional plotting 


“orm LAC 20} ae 5 


. 3's 


it can be done by reading the instructions for DRAW3D [2]. This report is available ' 
from computer operations. Also both NDT and NDP must be greater than zero to get a 
three-dimensional plot and if you try to plot a function with all zeros you will get 


a blank frame. 


The thermal history terms, which are a huge addition to the code, are described i 


below. 


The strain in the shell is composed of two components, one the kinematic (strains i 
due to shell deformation) and two the thermal (fictitious strains due to thermal ex- 


pansion or contraction). The kinematic strain is denoted as 


where i takes on the values $, @, $@- These are equations (2.1-2) in reference 


| 

i 

{1]. The thermal strain is denoted as | 
epi = CT - T) | 

| 

| 


where i takes on the values h) and 9, a is the coefficient of thermal expansion, 

tT is the reference temperature and T is the current temperature. It should be 

noted that £T is a function of $,0,z and t (the shell coordinates and time) so 

that, ‘ 


T = 1(f,0,z,%). 4 


The total strain becomes for $,6 and $@ coordinates 


eg + 2Ky + ag (T =") 


“$ 


e 
8 


+ - 
e.* 2K, Oy (rT ‘T) 


&h9 = &¢q * "Ko 


In terms of stresses we have 


Ey : E E E 
; -8 .( 4 8 : 
°% = - & + Yo Z €x 9 e °% + Vp - a.) (@ T,) 


where y=l - Ug as (ote that these are equations 2.1-3) of reference [1] and 


equations (2.1-3a&b) are incorrect in reference [1].) 


Briefly and simply, if the temperature increases with time we have positive 
thermal strain which produces a negative stress (compression). When the equilibrium 
equations are solved the compressive thermal stresses must be relieved by enema 
motion, hence the shell vibrates about some new equilibrium position that corres- 
ponds to a balancing of the stresses in the equilibrium equations. The code output 
for strain includes only the kinematic portion (i.e., what would be measured by a 
strain gage ) and the output for stress includes both the kinematic and thermal 


stresses as these are the actual stresses in the shell. 


The computational sequence is identical with that deveiorxed for the beam equa- 
tions used in TBYAM and TBEAM2 [3]. First in the input data one can specify all 
the shell material properties as a function of temperature through a tabular input. 
The data can be specified for eight different temperatures. Additionally the tem- 
peratures as a function of §,@ and z can be specified for eight different times, 


t, through the standard matrix read. 


During the response calculations the temperature for the time and spatial lo- 
cation are determined by a table lookup procedure using linear interpolation between 
the input data points. From this temperature value the shell material properties 
are then determined by table lookup. At this point the thermal terms in the strains 
and stresses, as outlined above, can be computed and then the response calculations 


proceed with the additional thermal terms. 


The temperature matrices are stored on a disc and only the values currently needed 
are pulled into core. The material properties array is stored in core so only a few 


additional storage locations are required for the thermal excitation terms. 
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Input Data Instruction Changes 


Cards 1 and 2 


a) As before two heading cards are read. 
b) As before the first card appears on the time-history plots. 
c) Now the first 48 spaces on the second card appear on the three-dimensional plots. 


a) As before both cards appear on the printed output. 


The meaning and use of I'EMP is different, now: 
a) ITEMP = O no thermal terms considered 


b) O< ITEMP <8 then IYEMP indicates the number of times at whicn the temperature 
matrices are going to be input (i.e., if your temperature-history is describable 


by four time points, ITEMP = 4). 


Gard 16 (10.3, 15,. 5X; E1013, 15) 
Where now PLT, NPLT, SPLT, NPSP are read. 


PLT - plotting time intervals ~ secs. 


NPLT 


Q no time-history plots 
5 ininiature plots only 
7 full and miniature 

8 full only 


spatial plotting intervals w~ secs. NOT end NDP 
must both be 
greater than 
> O spatial plots zero. 


SPLT 
NPSP - O no spatial plots 


Cards 23 
If ITEMP = O no cards are read 
If ITEMP > O the following cards are read in &E10.3 format: 


a) TEMB (8) - temperature values which correspond to the material properties 


below; the temperature values must be in increasing steps. 


Jf 


ge 
b) NUPD(8) - values of NUP corresponding to above temperatures 
c) wurp(8) - values of NUT corresponding to above temperatures 
ay Exp(8) - values of EX corresponding to above temperatures 
e) #rp(8) - values of ET corresponding to above temperatures 
bie GGD(8) - values of GG corresronding to above temperatures 
&) $1p(8) - values of Sl corresponding tc above temperatures 
h) s2D(8) - values of Sf corresponding te above temperatures 
i) slep(8) - values of S12 corresponding to above temperatures 


3) ALPHXD(8) 
k) ALPHTD(8) 


values of ALPHAX corresponding to above temperatures 


values of ALPHAT correspondi:.g to above temperatures 


ihe SzpD(8) - values of SZ corresponding to above temperatures 

m) HD(8) - values of H corresponding to above temperatures 

n) §iCD(8) - values of SIC corresponding to above temperatures 
eo) siTp(8) - values of SIT corresponding to above temperatures 
p) secp(8) - values of SeC corresponding to above temperatures 
qa) SeTD(8) - values of S2T corresponding to above temperatures 
r) $12sp(8) - values of $128 corresponding to above temperatures 


If NSL = 2 (i.e., a two layered shell cards aa-qq are also read, for the outer shell 


layer). 

aa) NUTAD(8) - values of NUTA corresponding to above temperatures 
bb) NUPAD(8) - values of NUPA corresponding to above temperatures 
ce) EXAD(8) - values of EXA corresponding to above temperatures 

dd) ETAD(8) - values of ETA corresponding to above temperatures 

ee) GGAD(8) - values of GGA corresponding to atove temperatures 

ff) SiAD(8) - values of S1A corresponding to above temperatures 

ge) Sean(8) - values of S2A corresponding to above temperatures 

hh) 312AD(8) - values of S12A corresponding to above temperatures 
ii) ALFXAD(8) - values of ALPHAXA corresponding to above temperatures 
jj) ALFPAD(8) - values of ALFHATA corresponding to above temperatures 
kk) SZAD(8) - values of SZA corresponding to above temperatures 

11) HAD(8) - values of HA corresponding to above temperatures 

mm) S1ACD(8) - values of SLAC corresponding to above temperatures 
nn) SiATD(8) - values of SiAT corresponding to above temperatures 
00) S2acD(8) - values of SZAC corresponding to above temperatures 


SE 


* SuGvs 
pp) S2ATD(8) - values of SeAT corresponding to above temperatures 
qa) S12asp(8) - values of S12AS corresponding to above temperatures 


End of cards for the outer shell layer 


s) TIMT(8) - time values which correspond to tne temperature arrays below. 
The time values must be in increasing steps. A zero time and temperatvre 
array must be included. If you want a step change in temperature, zero 
time must correspond to temperatures before heating and then at time - Df 
“nput the tempera!nres atter heating. The number of values of TIM? input 


must equal ITSMP. 
t) VOXEL (TL aJeb)  metris format 
The above temperature array is read for each time input in TIM. 


I = 1, NDP+1+NFREE , ( coordinate 
J = 1, NDT+1 , @ coordinate 


i 


L = 1, NLAY3+1 (L = 1 is the inner fiber), thickness 


Cards 30 (1215, NUSP - Spatial Plot Array 
It NPSP > © two new cards (right after cards 29) are read. 


Where if NLSP(I) = 0 no plot is made 
or NLSF(I) = 1 a plot is made 


where ‘he NLSP array corresponds to 


NLSP(1) - radial displacement” 

NLSP(2) - meridional displacem 

NLSP(3) - circumferential displacement 

NLSP(4) - outer fiber merd. strain 

NLSP(5) - inner fiber merd. strain 

NLSP(6) - outer fiber circ. strain 

NLSP(7) - immer fiber circ. strain 

NuSP(8) - outer fiber shear strain : 
NLSP(9) - inner fiber shear strain 

NLSP(10) - outer fiber merd. stress 


NLSP(11) inner fiber merd. stress 


x s : 
For the radial displacement you will get 3 plots (includes max and min cross-sections); 
all others are just 1 plot. 


9 


NLSP(12) 
NLSP(13) 
! NLSP(14) 
NLSP(15) 
NLSP(16) 
NLSP(17) 
NLSP(18) 
NLSP(19) 
NLSP(20) 
NLSP(21) 
NLSP(22) 
NLSP(23) 
NLSP(24+) 


’ 


= * 


outer fiber circ. stress 
inner fiber circ. stress 
outer fiber shear stress 
inner fiber shear stress 
merd. moment resultant 
circ. moment resultant 
shear moment resultant 
merd. stress resultant 
circ. stress resultant 
shear stress resultant 
radial pressure 

merd. pressure 


circ. pressure 


Don't forget the @FIN card at the end. 


we we 
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susJect: SUPPLEMENT II - USER'S GUIDE TO THE SHORE CODE {1} 


This supplement is tc provide formal documentation of the restart capability that has 
been added to the SFORE code. The production version of the SHORE code is available on 
the file PGU*SHOREP. 


The main new feature is the restart capability. Other changes include details of the 
plotting for better use of storage, that should not really affect the user and slight 
changes in the ronlinear formulation of the equilibrium equations. This later item 
should only produce very slight changes in answers; if you observe major changes in 
answers please notify me. The date of the last revision is now printed out on the first 


page of the output. If the code revision date changes it would be advisable to obtain 
a new listing. 


The input data is identical to that described in {1] and [2] except for two changes in 
Card 15 which are necessary if you want to restart. 


Card 15 (5E10.3, I5) DT, TT, PRINT, ESPT, START, NRST 


DT - integration time step (sec) 
TT - total time (the time span for this run) (sec) 
PRINT - printing time intervals (sec) 
ESPT - time left in run before plotting will start 
START - the time at which this run begins (sec) 
WRST - restart flag = -1 starts from a restart tape 
O nothing 
l writes a restart tape 
2 starts from a restart tape and 
writes a restart tape 


If you are going to write a restart tape it is wise to include the two following @MSG 
cards just before @xoT. 


€MSG N@TE TAPE T¢ BE ASSIGNED DURING RUN 
@MSG,I N@TE TAPE Té BE ASSIGNED DURTNG RUN 


If you are starting from a restart 
A 


you must read in the input data deck with the 
appropriate chanses in Card 15. e 


following control cards are needed, 


- 
Od 


@GASG,T RESTART., 16N, Xxx 
GUSE 16, RESTART, 


where XXXXX is the tape number oi the restart tape previously generated. 


NOTE: All tapes are 9 track 
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From one run to another you can change the following input data: 


1) Anything on CARDS 1 & 2 (heading) 

2) Anything on CARD 5 (print out flags) 

3) IMD on CARD 6 

4) Anything on CARD 9 (fracture stresses) 

5) Anything on CARD 12 (fracture stresses) 

6) Anything on CARD 15 except DT (DO NOT CHANGE DT)! 
7) Anything on CARD 16 (plotting times and flags) 

8) Anything on CARD 22 pressure data* 

9) Anything on CARDS 23 cemperatuve data* 
10) Any of the plotting data on CARDS 28, 29 or 30. 


*Avoid a discontinuity between runs, but you can add information for greater 
times that may occur after a restart. 


Finally, the restart tape will only be written if the job is successfully completed up 
to the start of the plotting operations, i.e. all the shell response calculations are 
completed. 

REFERENCES: 


(1) User's Guide to the SHORE Code, by Philip Underwood, IMSC-D244589, 30 Jan 1973 
[2] Supplement I - User's Guide to the SHORE Code, 19 Feb 1973 
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suBJECI: SUPPLEMENT III - USER'S GUIDE TO THE SHORE CODE [1] (A Status Report ) 


Introduction 


This supplement is to serve many purposes: 1) provide documentation on plane 
strain and stress ring analyses with the SHORE code and improved estimates of time 
steps, 2) a status report on what areas of the code have been studied this year, 

3) a report on areas that still need work and 4) a discussion keeping the code up 

to date and a better system for eliminating errors in the code. opefully the items 
discussed under 2 and 3 will spark some interest in users as the developer would like 
to farm out some development work in hopes that some fresh idcas will help some old 
problems. 


Discussion 


Ring Analyses 


If the two-dimensional version of SHORE is used for ring analyses (i.e. NDT=0 
and IAX=-1) one may set the input data up to do either a plane strain or plane stress 
analysis. The plane strain analysis would be representative of the midpoint behavior 
of a long shell loaded uniformly along the length. The plane stress analysis is repre- 
sentative of what is physically thought of as a ring (i.e. a very short free ended 
shell). It should be noted that the ring version of SHORE, SHORER, that is occasionally 
used is a plane stress analysis. The following input procedures should be followed for 
plane strain or stress. 


A) Plane Strain 


GG (shear modulus) = 0.0 
$12 (shear yield) >>> Sl & S2 


B) Plane Stress 


GG (shear modulus) = 0.0 

EX (¢-direction modulus) = 0.0 

NUP (¢-direction Poisson's ratio) = 0.0 
$l (G-direction yield stress) >>> S2 
$12 (shear yield) >>> S2 


Also to obtain the correct meridional strains for the plane stress case, some changes 
in subroutine STRESS are necessary. 


1) The statement DEPP = DEPPl + Z2*(... ), which is located two statements after 


357 CONTINUE, should be deleted and right after DEPT is calculated the statement 
DEPP = -NUTB*DEPT should be added. 
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2) Right after CALL INELST add 


SPT = 0.0 
sts = 1.0 
STS = SIGN(STS,DEPT) 


IF (L.EQ.1) EPP(II,JJ,1) + -NUTB*EPT(II,JJ,1) 
-(,.5-NUTB)*(EPT(II,JJ,1) - STS*S2B*RETB) 

IF (L.EQ.NLP) EPP(II.JJ.2 = -NUTB*EPT(II,JJ,2) 
-(.5-NUTB)*(EPT(II.57,2) - STS*S2B*RETP ) 


(If you are not interested in meridional strains there is no need to make the changes 
outlined in 1 and 2 above). 


Time Step Calculations 


In reference 2 more accurate expressions for computing the time step for stabile 
time integration are given. These expressions will eventually be placed in the SHORE 
code but until they are the user may want to check the time selected against the better 
expressions that are given below. 


C , Ss lia 
a= l(a) + (za) free 
and 


-1 
43 1 1 | 
a, 2S + 
fr he age (rap)* 


where At. is based on extensional frequencies and At, is based on flexural frequencies; 
at, is usually smailer than At, and the smaller At must be used. 


Also; h = thickness 

12 2 

C* = Eg/p(1-v*) 

and C,°= Eg/p(1+y) 
Status 


Status may not be the correct word for this section but somehow I wish to convey 
+o the user a couple of things that are considered to be deficiencies in the code and 
several areas of study that would add to the usefulness of the code are discussed. 


Sle 
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Deficiencies 


Two areas that have not been treated adequately are the spatial grid and the ( 
boundary conditions for nonlinear problems. 


The spatial grid problem pertains to the fact that all the stresses are needed 
at the same pivotal point to solve the plasticity problem. As pointed out in the i 
manual the grid is just shifted so all the stresses are computed at the same pivotal 
point. This shift reduccs the accuracy of the shear expressions. Indeed if one put 
a torsional load on the end of the shell the SHORE code would not obtain a very correct 
answer; the result would resemble the case of an axial end load shown in reference 3. } 
By averaging the derivatives in a certain manner this problem can be minimized to 
some extent; this has been done in the present SHORE and hopefully even better in a 
newer version of t1e code now being checked out. As a sidelight the nonlinear kinematic \ 
terms are also affected by the manner in which the derivatives are obtained. During the 
past year the nonlinear portion of the code was completely rewritten and to the best of 
my knowledge at this point in time is as accurate as can be considering the grid and 
shell equations used in the code. One interesting approach to this problem somewhat 
follows the method used in multi-dimensional wave propagation (references 4 and 5). 
Here one would use u, v, w and w’ as independent variables with u, v, and w’ at the 
corners of a quadrilateral grid with w at the geometric center. This of course implies 
a rotary inertia term to obtain the w’ and this may involve smaller time steps than are 
now required. Anyway, the spatial grid problem has been briefly outlined, I am still 
working on it, and the user should be very careful of problems involving direct inplane 
shear such as applied torsion. 


Now, on to boundary conditions; as many of you are aware the boundary condition 
solution fails to converge for large nonlinearities. This same problem also occurred 
in the beam codes and was solved by including the nonlinear term in the initia. guess 
for the displacement at the boundary. When the same approach was tried in the SHORE 
code two problems occurred. One is that shells are complicated compared to beams and 
second, the beam code uses total displacements and the shell code uses displacement 
increments. The nonlinear problem is much easier to solve with total displacements 
so a new version of SHORE based on total displacements is being checked out. When it > 
is running satisfactorily the addition of the nonlinear terms into the boundary conditions 
will again be tried and hopefully the problem of boundary condition convergence will be 
cleared up. 


Neglected Study Areas 


Various areas that would add to the SHORE code are briefly discussed. Hopefully 
the user may have some interest in these areas and can contribute to the evolution of 
the code. 


Storage - A more efficient use of core would increase the size of problems that 
could be studied and probably reduce costs for small problems that do not require 
large core storage. 


Plasticity and Fracture - The plasticity theory used should be improved to more 


accurately represent actual stress-strain relationships. Preliminary studies along 
the lines of references 6 and 7 have proven interesting and someone may be interested 
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in continuing them. Reference 8 provides a significant step in evaluating engineering 
plasticity and it also provides many areas for continuing work. Work in fracture is 
wide open to anyone with ideas. 


Segmented Shells - Once a boundary condition technique that is "fool proof", at 
least to some extent, is available a natural extension is to segmented shells. This 


would enable the study of discontinuity stresses at the junction of various shell 
geometries. 


Time Step - A criteria for the critical time step for nonlinear equations is 


needed to eliminate the need for runningproblems with very small time steps just to 
ensure stability. 


Rings and Stringers - A provision for analyzing ring and stringer stiffened 
shells would be valuable. 


Dynamic Relaxation - This techniques provides a way to obtain static solutions 
from a code that uses an explicit central difference time integration scheme. Some 
ground work has been done but there are still many unanswered questions before this 
technique will be useful as a production analysis method. 


Code Updating 


In the future new versions of the code will be announced on the print out of the 
current production version of the code. The new version will be described and the 
file name given. The old and new version of the code will be available simultaneously 
for approximately one month. If for any reason the user wants the ole version it will 
be their responsible to make their own file for it. After a new version is announced 
I would appreciate it if users would start using it for non rush jobs (if there are such 
jobs at Lockheed). If no ditriculties are encountered start using the new version 
for everything. 


Over the years I have learned that it is impossible to write a large computer 
code without including errors. To help me eliminated errors as fast as possible call 
me if the error needs urgent repair. Otherwise send me the run with the error and I 

will try to repair the code. Actually in either case if you send me the run with the 
error via the computer courier it is a great help. Secondly, if th-re are parts of 
the code the user feels are useless let me know and if there are things that would 


be more useful let me know. Hopefully, I can increase the useful ané@ cut out the 
useless. 


>). 

Pet lL eoreed, 

Philip G. Underwood 

Structural Mechanics Laboratory 
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